function [t,p,w]=spherequad(nt,np)

%SPHEREQUAD  Generate Gauss quadrature nodes and weights for numerically
% computing spherical surface integrals.
%
% [T,P,W]=SPHEREQUAD(NT,NP) computes the product grid nodes in
% r, theta, and phi in spherical and the corresponding quadrature weights
% for a sphere of radius 1. NR is the number of radial nodes, NT is
% the number of theta angle nodes in [0,pi/2], and NP is the number of phi 
% angle nodes in [0, 2*pi]. 
%
%
% Written by: Greg von Winckel - 04/13/2006
% Originally written for volume integration.
% Contact: gregvw(at)math(dot)unm(dot)edu 
% URL: http://www.math.unm.edu/~gregvw
% Modified by Dimitrios Karampinos for surface integration.


[x,wt]=rquad(nt,0); 
t=acos(x); wt=wt*1;     % theta weights and nodes (mapped Legendre)
p=2*pi*(0:np-1)'/np;        % phi nodes (Gauss-Fourier)
wp=2*pi*ones(np,1)/np;      % phi weights
[tt,pp]=meshgrid(t,p); % Compute the product grid
t=tt(:); p=pp(:);

w=reshape(wp*reshape(wt,nt,1)',nt*np,1);

